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In this work, a typical thermoacoustically driven pulse tube cooler as a no-moving part device has been 
investigated by a numerical method. A standing wave thermoacoustic engine as a prime mover in cou¬ 
pled with an inertance tube pulse tube cryocooler has been modeled. Nonlinear equations of unsteady 
one-dimensional compressible flow have been solved by the finite volume method. The model presents 
an important step towards the development of nonlinear simulation tools for the high amplitude ther¬ 
moacoustic systems that are needed for practical use. The results of the computations show that the 
self-excited oscillations are well accompanied by the increasing of the pressure amplitude. The necessity 
of implementation of a nonlinear model to investigate such devices has been proven. The effect of APAT 
length as an amplifier coupler on the performance of the cooler has been investigated. Furthermore, by 
using Lagrangian approach, thermodynamic cycle of gas parcels has been attained. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Thermoacoustics have been meant as a combination of thermal 
effects and sound by Rott [1 . He developed the mathematics 
describing acoustic oscillations in a gas in a channel with an axial 
temperature gradient, with lateral channel dimensions of the order 
of the gas thermal penetration which are much shorter than the 
wavelength. According to the sound wave, thermoacoustic devices 
can be categorized as standing wave and traveling-wave systems. 
According to the energy conversion, there are two types of ther¬ 
moacoustic effects: one is the acoustic oscillation powered by heat 
energy in thermoacoustic engine (prime mover), and the other is 
the heat flow driven by acoustic power in thermoacoustic cooler. 
In thermoacoustic engines, thermal energy from an external source 
is converted into kinetic energy in the form of acoustic waves. 
Energy conversion occurs in the stack (or regenerator) which 
smoothly spans the temperature difference between the hot heat 
exchanger and the ambient heat exchanger. During the last decade, 
numerous attempts have been proposed to investigate and to opti¬ 
mize the thermoacoustic engines [2-7]. Parametric studies have 
been performed for investigation of influences of working liquid 
[8,9 , stack geometry and resonator length [10], resonator diameter 
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[11], and resonance tube geometry shape [12] on performance of 
thermoacoustic systems. Furthermore, some general optimizations 
have been made 13-16 . 

Linear thermoacoustics has been expounded and summarized 
systematically by Swift in 1988 [17 . The theory supposes that 
the time dependence of all fluctuating quantities, such as pressure, 
velocity, temperature and so on, are purely sinusoidal. It is unable 
to display the thermodynamic cycle and the mechanism of ther¬ 
moacoustic conversion clearly. Besides, the linear thermoacoustics 
is suitable only for stationary thermoacoustic system, and cannot 
reflect its varying process. It is incapable of dealing with the 
nonlinear thermoacoustic phenomena, such as amplitude satura¬ 
tion, and frequency shift, which have been observed in lots of 
experimental thermoacoustic apparatuses. 

Nonlinear thermoacoustics as a branch of nonlinear acoustics 
began to be developed in the last decade. Karpov and Prosperetti 
from John Hopkins University have done a lot of outstanding work 
in this field [18-21 . They carried the expansion to fourth order in 
the perturbation parameter, and obtained explicit results for the 
initial growth, nonlinear evolution, and final saturation. Then, the 
dependence of the saturation amplitude upon several important 
parameters was illustrated. They also investigated the model 
performance on several examples, including a prime mover, an 
externally driven thermoacoustic refrigerator, and a combined 
prime mover/refrigerator system. Both their model and the present 
model are one-dimensional nonlinear model. However, there are 
some differences. In the present work, a model based on the finite 
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Nomenclature 


A cross-sectional area of gas flow (m 2 ) 

A v heat transfer area per unit volume (m” 1 ) 

C F Forchheimer’s inertial coefficient 

C P isobaric specific heat (J kg -1 K” 1 ) 

C s specific heat of solid (J kg -1 K” 1 ) 

C v isochoric specific heat (J kg” 1 K” 1 ) 

d diameter (m) 

d h hydraulic diameter (m) 

d g plate spacing (m) 

d w wire diameter (m) 

F general function 

/ friction factor 

j imaginary unit,=(-l) 1/2 

I< Darcy permeability (m 2 ) 

k thermal conductivity (W m” 1 K” 1 ) 

L length (m) 

l mesh distance (m) 

m mass (kg) 

rh mass flow rate (kg s” 1 ) 

Nu Nusselt number 

n number of screens per inch (0.0254” 1 m” 1 ) 

Pr Prandtl number 

q ext external heat load per unit volume (W m” 3 ) 
R gas constant (J kg” 1 K” 1 ) 

Re Reynolds number 

T temperature (I<) 

t time (s) 


u velocity (m s” 1 ) 

V volume (m 3 ) 

x longitude coordinate (m) 

Z compressibility factor 

^{} real part of 

Greek letters 

a heat transfer coefficient (W m” 2 K” 1 ) 

P opening area ratio of screen 

ILi dynamic viscosity (kg m” 1 s” 1 ) 

p density (kg m” 3 ) 

<j) porosity 

co angular frequency (rad s” 1 ) 

Subscripts 

I node number 

i face number 

S solid 

Superscripts 

(dot) time derivative 

Abbreviations 

PTC pulse tube cryocooler 

SWE standing wave engine 

SS stainless steel 


volume method with upwind scheme is developed to solve ID 
equations of a system with non-perfect gas (via considering a com¬ 
pressibility factor) as the working fluid. In many combined systems 
such as thermoacoustically-driven pulse tube cooler, high ampli¬ 
tude oscillations occur so that high order harmonics must be con¬ 
sidered to model these systems. The present model is not limited 
by a specified order of approximation. 

Thermoacoustically-driven pulse tube cooler (TADPTC) has the 
advantage of structure simplicity with no moving mechanical 
components and thus promises a potentially long lifetime. High 
frequency operation of such a system can lead to a much-reduced 
size, which is very attractive in small-scale cryogenic applications. 
Godshalk et al. 22] established the first high frequency TADPTC in 
1996. However, the cooler just reached a lowest no-load tempera¬ 
ture of 147 I< with the frequency being 350 Hz. 

Before 2005, pulse tube cooler was directly coupled with the 
thermoacoustic engine and obtainable pressure ratio for the cooler 
was limited by the capability of the thermoacoustic engine. In 2005, 
Dai et al. [23] proposed the concept of acoustic amplifier, which is 
actually a long tube connecting the engine with the pulse tube 
cooler. They established and studied a 300 Hz TADPTC utilizing 
the invention of acoustic pressure amplifier tube (APAT) to couple 
the standing-wave engine (SWE) with the pulse tube cooler (PTC). 
APAT is simply a thin circular tube with a length smaller than 1 /4 
of wavelength to connect the cooler to the engine. Their theoretical 
calculation showed that suitable length and diameter of the tube 
can lead to a pressure wave amplification effect which means that 
the pressure wave amplitude coming from the thermoacoustic 
engine can be much amplified to drive the pulse tube cooler. 

With modifications, lowest no-load temperatures of 95 I< [24], 
79.6 K [25], and 77.8 [26] were sequentially reached. Later, a low¬ 
est no-load temperature of 69.5 I< has been reached as reported on 
the work of Zhu et al. [27 . Comparing with the SWE previously 
used, Zhu et al. made two changes. Firstly, the inner diameters of 


the hot heat exchanger, the stack and the ambient heat exchanger 
were increased to increase inlet volume flow of APAT and the 
amplification ratio. Secondly, a cone-shaped resonator substituted 
the isodiametric resonator to decrease the viscous resistance. Zhu 
et al. [28] have investigated the influences of different parameters 
such as the engine stack lengths, the geometry of the resonator, the 
length of the APAT as well as the mean pressure on the overall 
performance of the device, experimentally. Most recently, Wang 
et al. [29] obtained a cooling power of 1.04 W at 80 K and a no-load 
temperature of 63 I< with 500 W heating power. 

The present work reports on the investigations of a 300 Hz 
frequency thermoacoustically driven inertance tube pulse tube 
cryocooler by using a numerical nonlinear model. The developed 
numerical model is based on nonlinear one-dimensional unsteady 
compressible flow equations solved by the finite volume method. 
APAT length as a key design parameter and its effect on the pres¬ 
sure ratio, the hot temperature of SWE, the cooling temperature 
of PTC, mass flow rate, and pressure amplitude are investigated. 

2. Physical model 

2.1. Geometry of the device 

The geometry of the device including the heat exchangers, 
stack, resonator, and buffers for SWE, and including the heat 
exchangers, regenerator, pulse tube, straightener, and inertance 
tubes, for PTC, and also including APAT (Acoustic Pressure Ampli¬ 
fier Tube) as a coupler, is shown in Fig. 1. 

2.2. Governing equations 

Because the mass flow, momentum flow, and the enthalpy flow 
in the longitudinal direction is very greater than those of lateral 
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Warm Ambient 


Hot Heat Heat Ambient Buffer 



Fig. 1. Schematic of a standing wave thermoacoustic engine coupled with a pulse tube cooler. 


directions are, a one-dimensional model can describe the perfor¬ 
mance of the system, acceptably. Thermal non-equilibrium model 
for porous parts (stack, HXs, and regenerator) is used. Governing 
equations include continuity, momentum equation, energy equa¬ 
tion, and equation of state for gas, and energy equation for solids 
of porous parts as follows: 

Continuity equation: 

l(pt) + 2- {P u$) = o (i) 


Momentum equation: 

g(^u) +l x (P^) = -4>f x -^u + C ^u |u|) (2) 


Energy equation for gas: 


d 

Ft 


P(j) ( CyT + 




f 


purj) yCpT + 



d_ 

dx 



+ (XAy(T s — T) 



Because of the existence of a cryogenic temperature region in the 
pulse tube cooler, a compressibility factor is considered in the 
description of the gas state. Hence, equation of state for gas is: 


P = ZpRT 



Energy equation for the solids: 

f) d f dT \ 

^[(1 - $)p s CsT s ] = m ((1 - 4 >%+ «A V (T - T s ) + q ext (5) 


Last term of Eq. (5) is the external heat per unit volume that is 
loaded to the warm heat exchanger and the cold heat exchanger. 


2.3. Boundary and initial conditions 


2.4. Empirical relations 


It can be easily shown that the gas-solid heat transfer area per 
unit volume for each medium (tubular parts, parallel plate stack, 
and wire mesh screens) can be calculated as follows: 


A 


v 


4 (p 


( 6 ) 


Porosity for tubular parts is 1 and for other parts can be easily cal¬ 
culated by using geometrical relations. Hydraulic diameters are pre¬ 
sented in Table 1. 

The friction forces in the momentum equation are formulated 
based on porous medium characteristics and are generalized for 
all parts. For the parts that Darcy permeability and Forchheimer’s 
inertial coefficient have not been defined, friction factor is related 
by using the following relation: 


dp 

dx 




c F p 

Vk 




Table 1 presents empirical relations of friction factors for different 
geometries used in the construction of the device. 

We notice that inside SWE, minor losses must be considered 
carefully as well as the major friction. Nevertheless, inside the 
pulse tube cooler, mentioned frictions are small in comparison to 
regenerator Darcy and Forchheimer’s effects. 

The heat transfer coefficient is defined as follows: 


a = 


kNu 


d\ 


( 8 ) 


Nusselt numbers for different parts are tabulated in Table 1. Rey¬ 
nolds number is: 


Re = 


p\}¥h 


( 9 ) 


No-penetration (zero velocity) condition must be satisfied at 
the four end boundaries (two ends of SWE and two ends of PTC). 

Initial conditions include zero velocity, charge pressure, and the 
ambient temperature. Because of the large heat capacity of the 
solid matrix, it takes thermal behavior a long time to reach periodic 
(cyclic) steady state. Heat transfer area of the parallel plate stack is 
very smaller than that of the regenerator. Thus, the temperature of 
the parallel plates is the last parameter which will reach the ther¬ 
mal periodic steady state. Since only the periodic steady state solu¬ 
tion is sought, the initial conditions for temperature are given as 
close as possible to the final solution. 

The program developed in the present model can let the gas 
oscillate by itself in acoustic engine because of numerical noises 
produced in the calculation procedure. However, an initial small 
mass flow rate is given to the program to reach steady conditions 
sooner. 


3. Numerical simulation 

To solve the Eqs. (l)-(5), the computational domain is divided 
into several finite volumes, and then the equations are integrated 


Table 1 

Empirical relations used in the model. 


Part 

L (cm) 

d (mm) 

Solid medium 

Hot buffer 

5 

50 

Tube 

Stack 

6 

50 

SS parallel plates 

Resonator 

70 

29 

Tube 

Ambient buffer 

29 

59 

Tube 

APAT 

50 

4.3 

Copper tube 

Regenerator 

3 

10 

600# SS screen 

Pulse tube 

4 

6 

Copper tube 

Inertance tube 1 

30 

2 

Tube 

Inertance tube 2 

110 

4.3 

Tube 
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i+1 


i+2 





Mass and energy control volumes 

Momentum control volume 

Velocity and mass flow vector 

Node of temperature, pressure and density 


Fig. 2. Staggered grid used in the model. 


CPU time in grid independence study, the simulations were per¬ 
formed for only 100 cycles. Several mesh densities were tested 
so that for each part of the system, the number of nodes was 
increased and the results were checked. For the stack, the number 
of nodes is more than that of the other parts, because the large 
temperature gradient of the stack makes larger false numerical dif¬ 
fusions in coarser grid. This will be sensible, especially in the calcu¬ 
lation of work or energy flows. The cycle-averaged acoustic power 
at the engine outlet (APAT inlet) was calculated according to Eq. 
(15) for the 100th cycle. Fig. 3 shows the acoustic power for differ¬ 
ent number of grids. Note that the presented values of Fig. 3 are not 
steady values, because the flow field solution requires about 1000 
cycles to reach steady periodic state. The calculated steady acoustic 
power at the engine outlet is 61.8 W. With the total number of 
nodes of about 350, it can be inferred that the results are indepen¬ 
dent of the grid system. 


over each corresponding volume (J-cv for continuity, energy and 
state equations and i-cv for momentum equation). A staggered 
grid is used for spatial discretization (Fig. 2). Integration of conti¬ 
nuity equation over a control volume gives: 

^ + -rhi = 0 (10) 


d(mu)j 

dt 


+ m !+ i U/ + i - rtiiUi = 


-Aj (/>,(?/ — Pi- 1) 



A Xiifii 
i 



d_ 

Ft 


m 



+ m i+ 1 






+ 0CiAyjVi(T s j — Tj ) 



d_ 

Ft 


[msCsTs^ = 


LA 


S' l s 


djj\ 
dx 


), 


i+1 

Ts,l) + QexFl 



+ ttiAvjVi{Ti 



P I V l = Z l m l K[ l (14) 

The second order backward approximation is used for the time 
derivative except for first step calculations, that first order back¬ 
ward is used. The second order upwind scheme is used for the con¬ 
vective (the enthalpy flow) terms as well as for the momentum 
flow. The thermal conductive terms are approximated by central 
difference scheme. Equation of state is put into continuity equation, 
so a relation is approached for pressure as an explicit function of 
mass flow rates. 

Finally, the equations are changed into algebraic form. The solv¬ 
ing process is as follows: 


1. Solving the momentum equation. 

2. Calculation of the pressure, explicitly. 

3. Solving the energy equation of the gas. 

4. Solving the energy equation of the solid. 

5. Repeating steps 1 to 4 until one time layer is finished. 

6. Repeating steps 1 to 5 until the periodic steady state is 
achieved. 


4. Results and discussion 

4.1. Grid Independence , and validity of the model 

To analyze the grid independence of the solution, a series of 
simulation runs have been carried out. In order to reduce the 


W ac =f j)(PuA)dt (15) 

In order to validate the model and method of solution, a com¬ 
parison is made against the previously published experimental 
work of Zhu et al. [28]. The system, including stack 2, resonator, 
1, and APAT 1 of their work is modeled. The characteristics of the 
system have been listed in Table 2. Input heating power is 
1000 W that is loaded to the warm heat exchanger. 

The time step size of 5.0 x 10 -5 s (i.e. each cycle is divided into 
66 steps, approximately), and a total node of 344 are used. For the 
parts with large temperature gradient such as stack, regenerator, 
and pulse tube, finer grid is used, whereas long parts such as APAT, 
resonator, and inertance tubes are discretized with coarse grid. 
Simulation run is performed for about 1000 cycles, so a steady 
periodic condition is reached. At each time step, 20 inner iterations 
are used for proper coupling between the equations and for better 
convergence. Running one simulation case by a computer with 
2.4 GHz processor and 1.0 GB RAM until the steady state takes 
about 1.6 h. 

Fig. 4 shows the cooling temperature versus cooling power of 
present results and the results reported by Zhu et al. [28 . The rela¬ 
tion between cooling power and cooling temperature (so-called 
cooling curve) is often nearly linear and the slope of the cooling 
curve is almost equal for two results. The present cooling curve 
is almost identical to that of Zhu et al. but less than that. The sim¬ 
ulated steady cycle-averaged no-load cooling temperature is 
63.3 K, which is lower than the experimental result 72 K. It must 
be noticed that the thermal parameters, especially, the cooling 
power are the most complicated parameters that have large errors. 
For other parameters such as operating frequency, flow rates, and 



Total Number of the Computational Nodes 


Fig. 3. Grid independence of the acoustic pressure. 
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Table 2 

Characteristics of the modeled system. 


Medium 

dh 

Friction equation 


Heat transfer relation 

Tube 

d 

Laminar: 

f 64 

J — Re 

Nu = 0.023Re0.8Pr° 35 



Turbulent: 

f 0.3124 

J ~ Re 025 


Wire screen 

<t> rf . 

(1-0) u wire 

0 1 2 33.6n0 20 3 4 C f 0.337n0 2 

k ~2 ip yr< ~ ft 2 


Nu = ( 1 + 0.9 9(RePr) 0 - 66 )^- 76 [30] 

Parallel plate 

2d e 

Laminar: 

f _ 96 

J — Re 

Nu = 1.3030Re a3271 Pr 1/3 [31] 



Turbulent: 

-L = 2/og,„(0.64Rev(f)-0.8 




Fig. 4. Comparison of cooling temperature versus cooling power for present results 
and the results of Zhu et al. [11]. 

pressure amplitudes, as presented in this paper, differences 
between the calculated and the experimental results are smaller. 
The main reasons of modeling errors are as follows: 

1. Secondary flows inside pulse tube are the most important losses 
that result in a reduction of the net enthalpy flow from the cold 
end to the hot end of pulse tube. One-dimensional models can¬ 
not predict the secondary flows. 

2. The empirical relations for friction factors and heat transfer 
coefficients have relative errors of about 10-15%. 

3. Conduction heat transfer through the walls of regenerator and 
pulse tube from the hot ends to the cold end of them has not 
been considered in this study. 

4. Radiation heat transfer from the hot ends of pulse tube and the 
solid matrix of regenerator to the cold end of them is a compli¬ 
cated phenomenon that has been neglected in the present work. 

4.2. Nonlinear aspects 

4.1.1. Self-excitation process 

One of the most challenging purposes is to capture the evolu¬ 
tion of the initial instability to steady state which the instability 
saturates to a finite amplitude. This is a strongly nonlinear effect 
and cannot possibly be traced by the models based on the linear 
theory. After transitioning from the initial velocity disturbance, 
the simulation yielded strong sustained pressure oscillations. 
Fig. 5 presents the time history of acoustic pressure at the engine 
outlet. The pressure amplitude tends to 0.21 MPa as it is reported 
previously [28]. An important capability of the present one-dimen¬ 
sional nonlinear model is precise tracing of the self-excitation pro¬ 
cess which a linear model has not it. Most of the previously 
reported time histories of pressure oscillations are obtained by 
CFD commercial softwares. It must be noticed that a smaller initial 
mass flow disturbance causes a longer time to reach steady peri¬ 
odic pressure oscillations. It was found that the flow field (pressure 
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Fig. 5. The time history of the self-excited acoustic pressure oscillations. 

and velocity) reaches the steady periodic condition, very faster 
than the thermal field does. The frequency which is calculated by 
the present numerical model is 303 FIz. It has an acceptable agree¬ 
ment with the experimental result of 300 FIz [28 . 

4.1.2. Steady state nonlinearities 

After attaining the steady periodic state solution, we concen¬ 
trate on steady nonlinearities of the results. Note that the linear 
model is based on a sinusoidal approximation of all variables and 
it is not permitted to be used when the pressure amplitude is 
greater than 10% of the mean pressure. For a general function 
F(x,t) that represents the velocity, pressure, density, and tempera¬ 
ture, Eq. (16) has been used for curve-fitting of the present results 
as functions of time and longitudinal coordinate. In most of nonlin¬ 
ear works, this form (up to a specified order) of unknowns was 
incorporated into the governing equations to derive x-dependent 
equations and to solve them. 

F(x, t) = F 0 (x) + 9?{f! (x)^} + 9 i{F 2 (x)e 2iC0C } 

+ 9J{F 3 (x)e 3ja)t } + ... (16) 

The gas temperature at the middle of the stack along with its 
first and second order harmonic approximations have been illus¬ 
trated in Fig. 6 to investigate nonlinearity of the flow field inside 
the standing wave engine. It is obvious that even the first order 
approximation acceptably fits the calculated data. The distinction 
between the calculated results and the second order approxima¬ 
tion is not clear. Inside the whole engine, the ratio of the pressure 
amplitude to the mean pressure is less than 6%. We investigate the 
other results inside the engine and conclude that the linear theory 
can predict oscillations satisfactorily, except a narrow region con¬ 
taining the minimum pressure amplitude point. As it is presented 
and discussed at the end of Section 4.3, some nonlinearities and 
jump to the second harmonic occur in this region of the resonator. 
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Fig. 6. The gas temperature at the middle of the stack versus time and its harmonic 
approximations. 

Inside APAT and the cooler, high amplitude oscillations occur. 
The results indicate that inside the cooler, the oscillations deviate 
from linear theory, because, the pressure amplitude is amplified 
by APAT so that it exceeds from 10% of mean pressure. The most 
important sources of nonlinearities are friction forces, including 
minor losses and Forchheimer’s inertia, and advective acceleration 
(in sections with very large diameter changes). Here, we present 
only two selected parameters: the mass flow rate at the cooler inlet 
and the gas temperature at the cold end of pulse tube. 

Fig. 7 shows the mass flow rate at the cooler inlet along with its 
first, second and third order approximations. First and second 
order approximations can not predict the mass flow rate oscilla¬ 
tions precisely. It can be inferred from Fig. 7 that the mass flow rate 
oscillates by third order harmonic. It is noteworthy that based on 
the conservation of mass, the cycle-averaged mass flow rate at 
any cross section must be zero, but due to the numerical errors, 
this criterion has never been satisfied and it tends to a very little 
steady value. 

The deviations are more obvious when the temperature field 
inside the cooler is investigated. The calculated results and the 
approximations (up to 5th order) of the gas temperature at the cold 
end are illustrated in Fig. 8. For better distinction, only first, third, 
and fifth order approximations are presented. The results show 



Fig. 7. The mass flow rate at the cooler inlet versus time and its harmonic 
approximations. 



Fig. 8. The gas temperature at the cold end of pulse tube versus time and its 
harmonic approximations. 

that a high order nonlinear thermoacoustic is necessary to model 
the temperature field of the gas inside the cooler. 

Hence, inside pulse tube, because of the high-pressure ampli¬ 
tude amplified by APAT, oscillations of the flow and thermal field 
are not sinusoidal, especially the temperature variations, so that 
implementing a nonlinear model is essential. 

4.3. Steady periodic state 

After the steady periodic condition is reached, we will concen¬ 
trate on one-cycle oscillations. To illustrate the order of magnitude 
of the gas that flows to the engine and to the cooler, Fig. 9 shows 
mass flow rates at branch faces at the moment of mean pressure. 
The magnitudes are close to maximum magnitudes of cyclic mass 
flow rates. Inequality of mass flow rates into branch control vol¬ 
ume causes a change in density of gas. 

Parametric study of APAT gives useful knowledge to optimize 
and investigate the flow phenomena inside the system. Inertial 
effects, pressure drop, and mass storage are the dimension depen¬ 
dent parameters that affect the amplitude and phase of pressure 
and mass flow rate, and the performance of pulse tube cryocooler 
subsequently. In order to study the influence of the APAT on the 
parameters, simulations were done by using an APAT with inner 
diameter of 4.3 mm and several different lengths. Curve fitting of 
the mass flow rate, velocity, temperature, and pressure has indi¬ 
cated that the oscillations are about sinusoidal inside the engine, 
but inside the cooler, deviations from sinusoidal behavior has been 
seen, especially in the temperature oscillations. 

APAT changes amplitude and phase of oscillating parameters 
while gas is flowing through it. Simultaneously, it decreases acous¬ 
tic power a little. Variations of mass flow rate at the engine outlet 
and the cooler inlet versus time for different APAT lengths are 
shown in Fig. 10. As APAT length increases, amplitude of mass flow 
rate increases. The ratio of mass flow rate at the engine outlet to 
mass flow rate at the cooler inlet increases as APAT length 
increases. APAT induces a phase lag nearly proportional to its 


j Density decreases \ 

105 g/s from AHX I ; | f 

l 4 g/s 1 

from APAT 

Fig. 9. Mass flow rates in branch at the moment of mean pressure. 
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Fig. 10. Mass flow rate versus time for different APAT lengths. 


length. The phase difference between mass flow rates of the engine 
outlet and cooler outlet is 7.6 deg for L apat = 0A m, 19.7 deg for 
L apat = 0.3 m, and 44.0 deg for L APAT = 0.5 m. 

Fig. 11 shows pressure oscillations for different APAT lengths. It 
is clear that a suitable APAT length can really lead to a more 
intense acoustic pressure wave available to the cooler than that 
provided by the engine. Increasing of APAT length has not signifi¬ 
cant effect on the amplitude of the pressure at the engine outlet 
so its value remains constant at the value of 0.21 MPa, approxi¬ 
mately. Nevertheless, the pressure amplitude at the cooler inlet 
and phase difference of the pressure wave at two ends of APAT, 
increase as APAT length increases. The phase difference between 
the pressure at the engine outlet and pressure at the cooler inlet 
for APAT lengths of 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6 is 32.70 deg, 
27.80 deg, 25.10 deg, 18.50 deg, 15.30 deg, and 0.90 deg, respec¬ 
tively. For better distinction, not all graphs of different simulated 
APAT lengths have been depicted in Fig. 10 and Fig. 11. 

For APAT lengths of 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6, phase differ¬ 
ences between the mass flow rate and pressure at the cooler inlet 
are 65.5°, 58.9°, 53.5°, 50.2°, 50.2°, and 50.2°, respectively, Whereas 
it is 27.3 deg for a typical compressor driven pulse tube cooler. 

APAT dimensions affect the flow field of the engine domain as 
well as flow field inside itself and the cooler. As it is illustrated 
in Fig. 12, volume flow rate inside the engine is influenced by APAT 



Fig. 11. Pressure oscillations for different APAT lengths. 


length. The volume flow rate inside the engine has a maximum 
value at the position of the resonator where the minimum pressure 
amplitude occurs, at x = 0.2 m from left end of the resonator. Inside 
the ambient buffer, the volume flow rate tends to zero linearly by 
steep slopes. Increasing of APAT length causes a larger amount of 
gas flows into PTC, thus the volume flow rate along SWE decreases 
a little. 

Three phenomena dominate the pressure field: 

(1) Amplification (inside APAT): The pressure amplitude 
increases. 

(2) Resonating (inside the resonator and inertance tube): a 
phase shift of pressure oscillation is produced to make the 
proper volume flow rate. 

(3) Friction (All parts; especially in porous media): The pressure 
amplitude decreases especially inside the regenerator. 

The distribution of the maximum and minimum of the pressure 
oscillation along the device are illustrated in Fig. 13. Inside tube 
parts such as hot buffer and pulse tube, there is a little change of 
magnitude of pressure amplitude. Inside APAT, the pressure ampli¬ 
tude increases and the velocity amplitude decreases. We can divide 
the resonator (also inertance tube) into three regions: left that 
includes about two third of the resonator length, mid-region that 
is a narrow region about a few centimeters, and right that about 
one third of the resonator length. In left region, the pressure ampli¬ 
tude decreases from the left (inlet) to the lowest pressure ampli¬ 
tude point, while the pressure phase does not vary so much. In 
the mid-region where the lowest pressure amplitude point is 
located, the pressure phase changes about 180° (see Fig. 14). This 
phase change is not provided by lagging but rather by jump from 
base harmonics to second harmonics. At x/L = 0.7 frequency is 
twice as large as the base frequency. The change of harmonics in 
Fig. 14 is evident. 

From the minimum pressure amplitude point to the right end 
(outlet), the pressure amplitude increases, while the pressure 
phase does not vary so much. 

4.4. Total performance 

From the previous experiences, a great performance could be 
easily obtained by using an APAT. Calculations indicate that a 
longer APAT length gives a better acoustic performance, but 
because of the safety requirements (the highest hot temperature 
of 923 I< in the present case), it is limited to a maximum value. 



Fig. 12. Amplitude of volume flow rate along SWE for different APAT lengths. 
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Fig. 13. The distribution of pressure oscillation limits along the device. 



Fig. 14. Pressure oscillations during one cycle at different (normalized) locations 
inside the resonator. 

Hot heat exchanger temperature for different APAT lengths is 
calculated and shown in Fig. 15. The experimental results are also 
presented for comparison. Obviously, when APAT length exceeds 
600 mm, highest temperature increases over 923 K. The difference 
between the numerical and reference values are less than 20 K, 
hence, an acceptable agreements are reached. 


The pressure ratio is defined as the ratio of the maximum pres¬ 
sure to the minimum pressure at the cooler inlet. The no-load cool¬ 
ing temperature and the pressure ratio are shown in Fig. 16. It 
must be mentioned that the graph of the pressure ratio is in good 
agreement with the experimental results except for small APAT 
lengths, because it must tend to 1.0 in the case that no APAT is 
used. It is obvious that there is an inverse relation between no-load 
cooling temperature and the pressure ratio. 

4.5. An extra challenge: Lagrangian approach 

By using the values of last calculated cycle in Eulerian coordi¬ 
nate, an additional Lagrangian approach is used to track the move¬ 
ment of different gas parcels to get their pressure, temperature, 
and velocity in a thermodynamic cycle. For the first time level, 
we fix Lagrangian coordinate (location of the particles) to Eulerian 
coordinate. Thus, all variables are same in the two coordinates. The 
gas parcel will move to a new position at the new time level as it 
oscillates. Assuming that the velocity of the gas parcel is constant 
during a time step, the new location of the gas parcel can be easily 
calculated. With the assumption that values are linear between 
two adjacent nodes, if the new position is not located in a grid 
point, pressure, temperature, and velocity of the gas parcel are cal¬ 
culated by linear interpolation of those parameters at two adjacent 
grids. After 66 time steps, the gas parcel will return to its original 
position at the beginning of the cycle. Therefore, thermodynamic 
cycles of the gas parcel can be obtained. Results show that different 
gas parcels undergo thermodynamic cycles at different tempera¬ 
tures. In this section, Lagrangian approach results of modeling of 
the device with specifications listed in Fable 2 are presented. 

For thermoacoustic engines, thermodynamic cycle of gas par¬ 
cels inside the stack is an elliptical path as seen from u-x diagram 
shown in Fig. 17 as well as T-x diagrams shown in Fig. 20. Fig. 17 
shows the velocity of three chosen gas parcels that are located in 
the stack. It is clear that the elements that are close to the cold 
end of the stack have more displacement and velocity variations 
than those are close to the hot end. The displacement of gas parcels 
is in order of 0.01 m. 

The acoustic power is proportional to the multiplication of pres¬ 
sure and velocity. APAT wastes the acoustic power due to wall fric¬ 
tion. The pressure increases along APAT, and velocity is decreased 
along it. It means that the gas parcels near the cooler inlet (APAT 
outlet) side have less velocity amplitude than those near the 
engine outlet as shown in Fig. 18. The gas parcels near the engine 
outlet, and near the cooler inlet have a velocity magnitude of 60 
m/s and 25 m/s. It implies that the linear thermoacoustic model 
cannot be used because of large inertial effects. The velocity varies 
linearly inside APAT. It can be seen from Fig. 17 that the gas parcels 
closer to the cooler inlet have less displacement. The maximum 
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Fig. 15. Hot heat exchanger temperature versus APAT length. 


Fig. 16. No-load cooling temperature & pressure ratio versus APAT length. 
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Fig. 17. Velocity of gas parcels inside the stack during one thermodynamic cycle. 


Fig. 20. Temperature of gas parcels inside the stack during one thermodynamic 
cycle. 



Fig. 18. Velocity of gas parcels inside APAT during one thermodynamic cycle. 



Fig. 19. Pressure of gas parcels inside APAT during one thermodynamic cycle. 


and minimum displacements of gas parcels inside APAT are 65 mm 
and 26 mm at the hot and cold ends of APAT, respectively. The 
gas parcel that is next to the cooler inlet, acts as piston of a 


reciprocating compressor for compressor-driven pulse tube cryo- 
cooler. Hence, its velocity is equivalent to the compressor fre¬ 
quency and its displacement is equivalent to stroke of the piston. 

Fig. 19 shows the pressure of three gas parcels inside APAT. The 
clockwise thermodynamic path means positive mechanical work 
(work input). The input power is equal to the area of the P-V cycle. 
Again, we can see that the gas parcels near the cooler inlet side 
have shorter displacement than those are near the engine outlet 
side of APAT. Along APAT, P-V diagrams become narrower, and 
the pressure amplitude becomes larger. It can be seen that the 
P-V diagram in the APAT inlet is tilted to the right largely, whereas 
gas parcels near the cooler inlet have a nearly symmetric path. 

Fundamentally, thermoacoustic devices can be either engines or 
refrigerators. It depends on direction of gas-wall heat transfer 
inside the stack in a standing wave device (or regenerator in a trav¬ 
eling wave device). For standing-wave devices, the direction of the 
heat transfer between the gas parcels and the solid surface near 
each end of the stack is decided by the critical temperature 
gradient. 

Fig. 20 shows the temperature of three local gas parcels as they 
oscillate along the stack. The cycle-averaged temperature of solid 
plates is also presented. For clarity, virtual linear temperature dis¬ 
tribution is shown. The cycle-averaged the temperature of the gas 
(in Eulerian coordinate) is equal to that of solid plates, approxi¬ 
mately. They can be estimated well with the following power 
function: 


T = T h 




Here, x is the longitudinal position from the hot end of the stack and 
L is the length of the stack. T c and T H are the cold end and the hot 
end temperatures of the stack, respectively. 

The clockwise direction of thermodynamic cycles is also illus¬ 
trated in Fig. 20. It can be clearly seen that the stack solid surface 
temperature is higher than the gas parcel near the hot end of ther¬ 
modynamic path, while the opposite is true near the cold end. The 
temperature of the parcel traces out a clockwise ellipse as a func¬ 
tion of time. 

As a gas parcel oscillates along the axis of the parallel plates, it 
experiences changes in temperature, caused by adiabatic compres¬ 
sion and expansion produced by the sound pressure and by 
gas-wall heat exchange. The gas parcels move to the hot side of 
stack while the pressure is rising. In this moment, the gas parcel 
absorbs heat from the adjacent plate wall. Then, it moves to the 
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cold side while the pressure is falling. Now, the gas parcel rejects 
heat to the plate. Imperfect gas-wall heat exchange is required 
in order to provide a time delay between the gas motion and gas 
thermal expansion or contraction. It results when the plate spacing 
is more than thermal penetration depth. 

5. Conclusions 

A nonlinear one-dimensional model for modeling of the flow 
and thermal field of a typical high frequency thermoacoustically 
driven pulse tube cryocooler has been developed. The finite vol¬ 
ume method has been used for discretization of governing equa¬ 
tions. The self-excitation process has been simulated very well. 
Results show that the present model has an acceptable accuracy 
in calculation of flow characteristics of the system. The effects of 
APAT on the performance of the device have been investigated. 
As an extra challenge, Lagrangian approach has been implemented 
to trace the thermodynamic cycle of the gas parcels. 
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